/**
 * @auth Ning
 * @purpose for Fourier Transform processing...
 * @date 2009-8-25
 * */	
package extendedmath.fourier
{
	import extendedmath.ComplexNumber;
	
	import utils.data.DataUtils;
		
	public final class FourierTransform
	{
		
		//START...
		
		// Discrete Fourier transform for one dimension
		public static function DFT(samples:Array):Array
		{
			
			var result:Array = new Array;
			
			var temp:ComplexNumber;
			
			for(var k:uint = 0; k < samples.length;k++)
			{
				
				temp = new ComplexNumber;
				
				for(var n:uint = 0; n < samples.length;n++)
				{
					
					if(k == 0)
					{
						
						trace(samples[n].multiply( (new ComplexNumber(-1,0)).pow(-2 * k * n / samples.length) ));
						
						
					}
					
					temp = temp.add(samples[n].multiply( (new ComplexNumber(-1,0)).pow(-2 * k * n / samples.length) ));
															
				}
				
				result[k] = temp;
				
			}			
							
			return result;
			
		}
		
		
		//IDFT...
		public static function IDFT(samples:Array):Array
		{
			
			var result:Array = new Array;
			
			var temp:ComplexNumber;
			
			for(var k:uint = 0; k < samples.length;k++)
			{
				
				temp = new ComplexNumber;
				
				for(var n:uint = 0; n < samples.length;n++)
				{
					
					temp = temp.add(samples[n].multiply( (new ComplexNumber(-1,0)).pow(2 * k * n / samples.length) ));
															
				}
				
				result[k] = (new ComplexNumber(1 / samples.length)).multiply(temp);
				
			}			
							
			return result;
						
			
		}
		
		//FFT...
		public static function FFT(samples:Array,retrorse:Boolean = false):Array
		{
			
			var result:Array = new Array;
			
			var store:Object = new Object;
												
			for(var k:uint = 0; k < samples.length;k++)
			{
											
				result[k] = processFFTElements(samples,k,retrorse,store);
				
				result[k] = retrorse ? (new ComplexNumber(1 / samples.length)).multiply(result[k]) : result[k];	
																		
			}
						
			return result;
			
		}
		
		//process FFT elements...
		private static function processFFTElements(samples:Array,k:uint,retrorse:Boolean,store:Object,parity:uint = 0):ComplexNumber
		{
									
			var result:ComplexNumber;
			
			var com1:ComplexNumber;
			
			var com2:ComplexNumber;
			
			var splitArr:Object;
			
			var wn:ComplexNumber;
						
			if(samples.length == 1)
			{
				
				result = samples[0] is Number ? new ComplexNumber(samples[0]) : samples[0];
												
			}else{
				
				var arrK:Object = store[[k < samples.length / 2 ? k + samples.length / 2 :  k - samples.length / 2,samples.length,parity]];
				
				if( arrK is Array)
				{
					
					com1 = arrK[0] as ComplexNumber;
					
					com2 = arrK[1] as ComplexNumber;
					
					result = com1.sub(com2);
					
					
				}else{
				
					wn = new ComplexNumber(-1,0).pow((retrorse ? 2 : -2) * k / samples.length);
									
					splitArr = DataUtils.splitArrayByParity(samples);
								
					com1 = processFFTElements(splitArr["even"],k,retrorse,store);
				
					com2 = processFFTElements(splitArr["odd"],k,retrorse,store,1);
				
					com2 = wn.multiply(com2);
				
					store[[k,samples.length,parity]] = [com1,com2];
				
					result = com1.add(com2);	
					
				}
								
			}
			
			return result;
			
		}
		
		
		//my fft...
		public static function myFFT(samples:Array,r:uint = 0):Array
		{
									
			var sampleLength:uint = samples.length;
									
			if(r == 0)
			{
				
				var z:uint = 1;
				
				while((z <<= 1) <= sampleLength)
				{
					
					r++;
					
				}
											
			}
			
			var sliceUnit:uint;
			
			var maxTempNum:Number;
			
			var tempArr:Array = samples;
			
			var storeArr:Array;
						
			// get coefficient array...
			
			var arrTemp:Array = new Array;
			
			//init arrTemp...
			for(var i:uint = 0;i < sampleLength;i++)
			{
								
				var k:uint = i;
				
				var c:uint = 0;
								
				for(var u:uint = 0; u < r ; u ++)
				{
					
					var d:uint = k % 2;
															
					c = c * 2 + d;
					
					k = d == 0 ? k / 2 : (k - 1) / 2;		
									
										
				}
												
				arrTemp[c] = new ComplexNumber(samples[i]);
								
			}
			
			var wn:ComplexNumber = new ComplexNumber;
			
			for(var w:uint = 0;w < r;w++)
			{
								
				storeArr = new Array;
				
				var length:uint = sampleLength >> (r - w - 1);
				
				var rb:uint = 1 << (r - w - 1);
								
				for(var b:uint = 0; b < rb;b++)
				{
				
					for(var a:uint = 0; a < length / 2;a++)
					{
						
						var index:uint = length * b + a;
						
						var odd:ComplexNumber = arrTemp[length * b  + length / 2 + a];
												
						wn.setTriangleForm(1,-2 * Math.PI * a / length);
						
						odd = wn.multiply(odd);
						
						var even:ComplexNumber = arrTemp[length * b + a];
						
						storeArr[index] = even.add(odd);
						
						storeArr[index + length / 2] = even.sub(odd);
						
					}
				
				}
								
				arrTemp = storeArr;
				
												
			}
			
				
			return storeArr;
			
		}
		
//		 public static function newFFT(TD:Array,r:int=0):Array {
//                //算法必需变量
//                var b,c,i,j,k,p:int;
//                var a:Number;
//                var WR:Array=new Array(c>>1);
//                var WI:Array=new Array(c>>1);
//                var FD:Array=new Array(c);
//                var XR1:Array;
//                var XR2:Array=new Array(c);
//                var XR:Array;
//                var XI1:Array=new Array(c);
//                var XI2:Array=new Array(c);
//                var XI:Array;
//                //为了减少AS3的计算量而增加的辅助变量
//                var t1,t2,t3,t4,t5,t6:int;
//                //计算傅立叶变换点数
//                if(r==0){
//                var z:int=1;
//                var l:int=TD.length;
//                while((z<<=1)<=l){
//                        r+=1;
//                }
//                }
//                c=1<<r;
//                //加权系数
//                for (i=0; i<(c>>1); i+=1) {
//                        a=-i*(c<<1)*Math.PI;
//                        WR[i]=Math.cos(a);
//                        WI[i]=Math.sin(a);
//                }
//                XR1=TD.slice(0,c);
//                for(i=0;i<c;i+=1){
//                        XI1[i]=0;
//                        XR2[i]=0;
//                        XI2[i]=0;
//                }
//                //蝶形运算开始
//                for (k=0; k<r; k+=1) {
//                        t1=1<<k;
//                        for (j=0; j<t1; j+=1) {
//                                b=1<<(r-k);
//                                for (i=0; i<(b>>1); i+=1) {
//                                        p=j*b;
//                                        t2=i+p;
//                                        t3=t2+(b>>1);
//                                        XR2[t2]=XR1[t2]+XR1[t3];
//                                        XI2[t2]=XI1[t2]+XI1[t3];
//                                        t4=XR1[t2]-XR1[t3];
//                                        t5=XI1[t2]-XI1[t3];
//                                        t6=i*t1;
//                                        XR2[t3]=t4*WR[t6]-t5*WI[t6];
//                                        XI2[t3]=t4*WI[t6]+t5*WR[t6];
//                                }
//                        }
//                        XR=XR1;
//                        XR1=XR2;
//                        XR2=XR;
//                        XI=XI1;
//                        XI1=XI2;
//                        XI2=XI;
//                }
//                //重新排序
//                for (j=0; j<c; j+=1) {
//                        p=0;
//                        for (i=0; i<r; i+=1) {
//                                if (j&(1<<i)) {
//                                        p+=1<<(r-i-1);
//                                }
//                        }
//                        FD[j]=Math.sqrt(XR1[p]*XR1[p]+XI1[p]*XI1[p]);
//                }
//                WR=XR1=XR2=XI1=XI2=null;
//                return FD;
//        }
//		
		
		//END
	}
}